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Abstract 

We consider a variational method to solve the optical flow problem with 
varying illumination. We apply an adaptive control of the regularization pa¬ 
rameter which allows us to preserve the edges and fine features of the computed 
flow. To reduce the complexity of the estimation for high resolution images 
and the time of computations, we implement a multi-level parallel approach 
based on the domain decomposition with the Schwarz overlapping method. 
The second level of parallelism uses the massively parallel solver MUMPS. We 
perform some numerical simulations to show the efficiency of our approach 
and to validate it on classical and real-world image sequences. 


Introduction 

In the last decades, the estimation of the optical flow has become a popular and 
central problem in computer vision m, m, m, ci)- It is involved, for example, 
in almost all the movies and pictures compression processes, or in the obstacle 
detection in the new smart cars and in robotics. The modelling and the detection 
of the motion in a scene involve numerous difficulties (e.g. the aperture problem, 
occlusions, ... m)- Among such difficulties, for the optical flow estimation, one 
with growing importance is the cost of the method in term of computation time 
(and the storage), which rises with the increasing resolution of the images due 
the technological devices progress. Up to now, there exist numerous methods for 
the optical flow estimation among which the Partial Differential Equations and 
particularly, the variational methods turn to be very efficient. They offer a complete 
framework which consists of mathematically founded continuous models, and a large 
number of numerical methods m, m, m, c])- They allow to cope with the ill- 
posedness of the optical flow problem, due to the aperture problem [7], by including 
a large range of regularization procedures. In this article, we consider a variational 
method based on the linear Horn and Schunck approach HO] but with a variable 
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regularization parauieter. This method introduced in [4] is proved to be an efficient 
approach to solve the optical flow problem in the sense that it is edge preserving 
and low cost (in term of degrees of freedom) [3]. 

A large part of the research works in the literature deals with the determination 
of optical flow under an assumption of constant illumination. This assumption is 
not satisfied in general and becomes a source of inaccurate estimation and serious 
limitations in many applications. Besides, modelling a varying illumination is quite 
complex and increases the ill-posed character of the problem. In this article, fol¬ 
lowing Gennert and Negahdaripour [8], we assume given an a priori law for the 
illumination variations and we introduce a supplementary variable, which models 
the illumination, in the system of equations. We show that the obtained partial 
differential equations system solved in the framework of the adaptive variational 
approach [4] is an efficient and innovative method for the optical flow estimation 
with varying illumination. 

A classical criticism against the variational methods is their ’’complexity" and 
their potential time consuming character, particularly when using unstructured 
meshes and finite element method discretization. The main contribution of this 
article is to propose and to validate an efficient massively parallel multi-level solver 
using domain decomposition method to solve the optical flow problem with varying 
illumination, following the variational adaptive approach proposed in [4]. We ob¬ 
tain the optical flow with an accurate estimation and a significant reduction of the 
time of computations. We validate our method on a classical benchmark and two 
real-world image sequences. 

In Section 1, we recall the Horn and Schunck model for the optical flow esti¬ 
mation and we give the system of equations to solve in the case where we allow 
illumination variations. In Section 2, we briefly recall the finite element method 
and we will rewrite the discrete optical flow system. We also define the adaptive 
strategy and give the resulting algorithm. In Section 3, we will consider the domain 
decomposition method with overlapping and the additive Schwarz method used to 
solve the algebraic problem. The Section 4 concerns the numerical simulations. We 
present some results obtained with our method, in particular, we give in detail the 
optical flow estimation for the so-called RubberWhale sequence to show the per¬ 
formances of the approach. We also perform the analysis of the computation time 
of the massively parallel algorithm. Finally, we will present two examples of the 
computation of the optical flow for real-world sequences. 


1 Optical flow problem 

We consider a sequence of two successive frames where fi c is the image domain. 
The intensity of a pixel {x, y) at an instant t is defined by the function 

I: If X [0, T] ^ R 

{x,y),t I{x,y,t). 

As it is usually done [2] , we use a convolution with a Gaussian kernel of standard 
deviation a to work with smoothed images. We define the smoothed image sequence 
by 

f{x,y,t) = {K„-kI){x,y,t). 
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The optical flow u = {ui^U 2 ) represents the vector field describing the motion of 
each pixel in the sequence. Following Horn and Schunck cni, we first consider the 
brightness constancy assumption. It means that the brightness intensity stays the 
same between two successive frames, 

= f{x + Ui,y + U2,t + l) ( 1 ) 


which is equivalent to say that the variation in time is null 


dt 


{x,y,t) = 0. 


Since the displacements are supposed to be small, we can assume that / is C^([0, oo[; M). 
So by using a first order Taylor expansion, we obtain the fundamental constraint of 
the optical flow 


dt dx dt dy dt 


( 2 ) 


The time derivatives of x and y represent the component ui and U 2 of the optical 
flow. Hence, with the notation /* = the equation becomes 


fxUl + fyU2 + /t — 0. (3) 

In this way, we have to determine two unknowns ui and U 2 with only one equation. 
This is the so-called aperture problem illustated in the figure In the example 1, 
in a local neighborhood we can only detect the vertical motion. In the example 2, 
it is only the horizontal motion that we can find. The example 3 is the only one 
where we can detect a diagonal motion. 



Figure 1: Representation of the aperture problem. 

To go through this ill-posedness, Lucas and Kanade [H] assume that the motion 
is constant in a neighborhood of size p. Contrary to this local assumption, Horn and 
Schunck propose cni a global approach to overcome the aperture problem. They 
introduce a regularization parameter a which acts as a penalizer and leads to a 
smoother flow field (the bigger a is, the smoother the flow is). The idea of Bruhn 
and Weickert [7) is, to combine both methods and minimize the functional 

/ Kp-k{fxUi + fyU2 + ftf + a{\Vuif+\Vu2f)dxdy (4) 

Jn 

where Kp is a Gaussian deviation of parameter p and a is a constant regularization 
parameter. 
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On real-world images, the assumption of the brightness constancy is no longer 
verified. Occlusions, shadows or glints don’t meet this constraint. Hence, the es¬ 
timation of the previous model is not accurate so we need to consider another 
assumption. Different approaches were proposed to model illumination variations, 
such as the assmption of the constancy of the gradient amplitude [6] or, in the case 
of color images, we can cite the work with variables less sensitive to such illumi¬ 
nation changes m- In this article, following Gennert and Negahdaripour [8], we 
consider a varying illumination obeying an affine transformation. This assumption, 
even covering a wide range of applications, might appear as a strong one. Moreover, 
it introduces a supplementary unknown, enforcing the ill-posedness of the optical 
flow estimation. To balance such potential shortcomings, we couple this modelling 
with an adaptive control of the regularization of the parameter associated to the 
new unknown. In this article, we restrict ourselves to smooth variations, keeping the 
parameter large enough. The assumption allowing a linear motion of the brightness 
intensity between the two images becomes 

f{x + ui,y + U2,t+l) = M{x,y,t)f{x,y,t) + T{x,y,t) 

where M is the multiplier and T is the translator. In our case, we suppose that the 
translator is negligible. The smaller the displacement is, the closer to 1 M is. In 
this way, we can set 


M = 1 + 5m (5) 

where 5m tends to zero when the displacement is very small. That gives the new 
equation of the optical flow 


fxUi + fyU2 + ft- fnit = 0 ( 6 ) 

where rtit is the derivative of M. For more details, we refer the reader to [H]. The 
optical flow problem allowing varying illumination consists finally of minimizing the 
functional 



* {fxUi + fyU 2 + ft- rntff' + a(|Vuip+|Vu 2 p) + X\Vmt\^dxdy. 


( 7 ) 


According to Euler-Lagrange equations, we have the system 


with 


r -div(AVU) + ApJJ = F in (8) 

I ^ = 0 on aff 
on 



a 


Ui 
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It is well known that taking a and A constant, even well chosen, leads to unde¬ 
sired oversmoothing and blurs the edges of an image. Thus, following the idea of 
( 0 , 0 ), we will consider the general setting where a and A are discontinuous and 
piecewise constant functions in order to prevent these smoothing effects. Indeed, 
as proved in [3], choosing a small value of a in regions where there are edges gives 
sharper edges, then an improved restitution of the motion and its segmentation. We 
refer to [3] for more details. In particular, we assume given fl = Ufli, a = 

i 

A = and the non-negative minimal value of a^. Then, we have the 

following theorem. 


Theorem 1.1 


The problem has a unique solution in 


The well-posedness is obtained from the Lax-Milgram lemma. By assuming that Tt 
is Lipschitz, we can work in the Hilbert space We define the norm 


1 

where \\Ap\J\\‘j^ 2 (^Q^= (ApU,U). Let denote Pm = max((aM 5 Am) and 

Pm = iiiin((a^, A^), a straightforward extension of ([3], proposition 2.7) gives the 

following proposition. 

Proposition 1.1 Let U he a solution of There exists C > 0 independent of a 
such that for Ua, the solution of the optical flow problem where the regularization 
is a piecewise constant function, the following inequalities hold 

\\UA\\p,A<C\\Alu\\mny (9) 

and 

1 

II u- Ua\\p,a< C ' \\a--VUA\\LHny (10) 


From now on, we denote U in place of Ua for brevity. 


2 Finite elements implementation 

2.1 Variational form 

In order to solve the system Q, we use the finite element method. To do so, we need 
to write this equation under its variational form. We multiply the first equation of 
(|^ by a test function (p, we integrate over Q and by using the Green formula on 
the first integral, we have 

[ AVWcpdyiT [ ApU(^dx= [ Fc^dx, \/cp e . (11) 

Jn Jn Jn 
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2.2 Discretization 

We consider a regular triangular mesh Th (figure . The estimation consists of 
solving a linear problem on each cell of the mesh. 



Figure 2: Representation of the initial mesh 7^. Each square represents a pixel of 
the image and consists of two triangles. 


We define the space of approximations by 

Vft = {V;, G Wh\KePi{Kf} 

where Pi{K) is the space of the linear functions on G Th- If we set Ap^h a finite 
element approximation of Ap^ the discrete problem of the optical fiow states 

{ Find U/i G V|, such that 

f AVVh-VVhdxdy+ [ VlAp,hUhdxdy = [ F^-Vhdxdy, VV;, G 

JTh JTh JTh 

( 12 ) 

We can show, by using the Lax-Milgram lemma that this weak formulation admits 
a unique solution. 


2.3 Adaptive regularization 


In this part, we are interested in controlling the parameter a. The regularization 
parameter is now considered as a piecewise constant function. This local choice of 
a is based on an a posteriori strategy analysis and was proposed by Belhachmi and 
Hecht [4]. The choice of the regularization is motivated by the fact that a small 
value is usefull to correctly approximate the Neumann boundary conditions on the 
edges of objects. However, it increases the maximum value of the optical fiow. So, 
in order to have a better estimation, we prefer a large regularization. The local 
choice of a allows to decrease its value on regions where we need a small a and keep 
a large value in the rest of the image. In our case, we have the additional unknown 
rrit with its regularization parameter A and in a first time, we keep this parameter 
constant. 

Since we want to locally choose the regularization parameter, we will have a large 
ratio so we will use the inequality (10) in the error indicator. The control of 


the regularization is done through an error indicator which is given for each element 

KeThhy 


Vk = + div(AKVUA,/i) 


1 


eesK 


Ae"/l|||[AVUA,/.' 


n. 
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where sk represents the set of all edges e of K. The diameter of K is noted Hk 
and the diameter of an edge e is hg. represents the normal vector from e, Ag is 
the maximum between the A of the two neighbors of an edge, and [.]g represents 
the jump over the edge e which means the difference between the outside and inside 
values. The error indicator r]x describes the finite element error and the model 
error. On the potential set of discontinuities, this value is large because is 

large. In fact, when the brightness is abruptly changing in an area, it means that 
we are close to an edge for the optical flow. So, to improve the solution, we decrease 
a from the two first components of A and the third component A stays constant. 
The decreasing formula for a is given by 



where k. is an arbitrary control parameter and ath is a threshold. In this way, if 
the relative error (measured with tjk) is greater than a given value rj we reduce 
in the value of a. On the other hand, if it is less than 77 the denominator is 
equal to one and so a remains unchanged. This local and adaptive control of a is 
implemented with the algorithm 

1 . Compute of a first approximation of the optical flow. This estimation is 
done on a cartesian grid where we have one cell per pixel. Define i = 0. 

2 . i = i 1. Build an adapted mesh with the metric error indicator. 

3. Update of ai{x,y) on T^. 

4. Go to step 2 . 

The convergence of this algorithm when the mesh size goes to zero is proved in [3]. 

3 Domain Decomposition 

3.1 Image decomposition 

In the case where we want to estimate the optical flow between two large images, 
we have implemented a domain decomposition. Indeed, as we have seen on the 
figure if we have large number of pixels on the picture, we have a large number 
of triangles in the mesh. In the finite element method, the resolution consists of 
solving a large system with a large number of degrees of freedom. The aim is that 
every Central Processing Unit (CPU) computes the optical flow of one part of the 
image (see figure]^. So, the number of linear systems to solve is reduced by a factor 
equivalent to the number of CPU used. The domain decomposition method allows 
us to obtain performances which are difficult to reach with classical variational 
methods and even worse with the adaptive process to control the parameter a. 
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CPU 2 


CPU 0 



CPU 3 


CPU 1 


Figure 3: Decomposition of an image for four CPU. 


Each CPU communicates the computed flow on the common boundary to his 
neighbor. So the more CPU we have, the more important the communication time 
is. More, the estimation of the optical flow is very sensitive at the boundary. For 
these reasons we want to optimize the number of pixels at the boundary of each 
part of the image. This is why we look for the largest ratio where area 

is the total number of pixels in the sub-domain and perimeter is the number of 
pixels on the boundary. In the figure we propose an example of a 48 x 48 grid 
and a decomposition for twelve CPU. We give the different values of this ratio with 
respect to every possible splittings. Because the grid is squared in this example, 
the ratio is the same for the symmetric decompositions. Finally, we can see that 
the ratio is better if we split the image in 4 x 3 (or 3x4). 




Figure 4: Evolution of the ratio on a 48 x 48 grid and twelve CPU. 

Since the optical flow estimation is sensitive at the boundaries, we use an additive 
Schwarz method to improve the estimation on the interfaces. This method requires 
an overlapping between the subdomains. 
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3.2 Model decomposition 

We note fti the part of the image corresponding to the CPUi and Ji the set of all 
indexes j which are neighbors to i. We define D Qj and Tij = dT^ij\dQj 

(See figure 1^. 



Figure 5: Example of notations for CPU i = 2. 

By using the additive Schwarz method to find an estimation of the optical 
flow in the part the problem related to @ is 


-div(A'=VU,^) + ApV'l = F in 


dn 

u? = on T 


0 on dQ.i\Ti^j, Vj e Ji 


(14) 




Vi e Ji. 


Thanks to this method the sequences (U^) converge to . The rate of con¬ 
vergence increases with respect to the size of Another convergent method 

presented by P.-L. Lions El exists with the Robin boundaries conditions. 


3.3 Multi-level parallel method 

To solve the system with the finite element method we use the sofware FreeFem+ + 
[S]. The default solver of this software for linear systems is UMFPACK (Unsym- 
metric MultiFrontal method) which allows to use non-symmetric matrices. The 
drawback of this library is that it can only solve problems smaller than 4GB. It 
means that we can’t treat parts larger than about 500 x 500 pixels per CPU so we 
need to cut the whole image enough in the case of large images. This implies to 
have a large number of CPU which is not always the case depending on the machine 
we work with. 

To overcome this issue, we use the MUMPS (MUltifrontal Massively Parallel 
sparse direct Solver) library which allows the resolution of sparse linear problems 
in parallel. The advantage of this solver is that it is not limited by the size of 
the problem, so we are able to use high definition sequences. It also ensures, to 
some extent, the scalability (see figure 13). However, it enforces to apply an LU 


decomposition to the matrix which can be long if this one is large (see figure 15). 
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Combining the domain decomposition method with the use of the MUMPS library 
(see figure [^, we can reduce the computation time of the high definition sequence 
estimation. 



Figure 6: Decomposition of an image for four groups of CPU. 

The image is split in the same way as above according to the maximal ratio 
between the perimeter and the area of the part. The number of the group in which 
belongs the CPU^ is given by 

group(CPUi) = i%(nbPart) 

where the binary operator a%b states for the remainder of the division of a by 6 and 
nbPart is the number of parts in the split image. The rest of the implementation is 
the same that the domain decomposition without MUMPS except that we act on 
a group instead of a single CPU. 

4 Numerical Results 

To test our algorithms, we use the RuhherWhale sequence (figure]^ given by the 
Middleburry website at www.vision.middleburry.edu/flow/, 



Figure 7: Frames 10 and 11 of the Rubber Whale sequence. 

Following their convention, we represent the estimated vector field thanks to 
a color map (figure which assigns a color to each vector with respect to its 
orientation and its norm. 
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Figure 8: Vector field and its corresponding color map. 


To validate our implementation, we compare our results with the exact solution 
also provided by Middelburry (figure [^. 



Figure 9: Ground truth solution of the Rubber Whale sequence. 


For optical flow problems, the accuracy of the method is usually evaluated by 
computing the Average Angular Error (AAE) given by 


AAE = acos 


_+ U2^hU2,e + 1 _ 

y i^l,h + '^2,h + l)(^l,e + ^2,e + 1) 


and the Endpoint Error (EE) given by 


EE = ^{Ui^h - + {U2,h - M2,e)^ 

where Uh = is an approximation of the vector field and Ue = ('i^i,e:'^2,e) 

represents the exact optical flow. Using the resolution of the system Q, our al¬ 
gorithm reaches an average angular error equals to 20.89 and an endpoint error 
equals to 0.38. In the litterature, the best angular errors go from 1 to 15 and 
the endpoint errors from 0.07 to 0.39 for the equivalent evaluation test case called 
Army (see the evaluation table at http://visioii.middlebury.edu/flow/eval/ 
results/results-_al .php). There are two reasons to explain that our error is 
large compared to these values. Eirst, we recall that EreeEem++ [9] uses unstruc¬ 
tured meshes and the computation of the angular error is not invariant with respect 
to the choice of the mesh. The other reason is that, even if we obtain a good 
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approximation of the vector field direction (see the results presented bellow), due 
to the large value of the regularization in the non-egde areas, we under-estimate 
the vector norms. There exists some iterative strategy to improve this value but 
since we are principally interested in improving the computational time and in the 
adaptation of the regularization parameter, we don’t use it. To validate the domain 
decomposition, we need to verify the convergence of the Schwarz method which 
means that 


I _ ijk+i I —0. 

k^oo 

We have split the image in four parts and tested the impact of the size of the overlap 
by using three different values, see figure 



Figure 10: Convergence of the Schwarz method for three different overlaps. 

We can see that the Schwarz method has a fast convergence and that the speed 
of convergence increases with the size of the overlap. However, if the overlap is 
large, the time of the construction of the matrix H, which is the longest part of the 
code, is large too. Hence, we need to choose the size of the overlap considering an 
optimal ratio between the rate of convergence and the speed of computation time. 
Therefore, we choose an overlap of five pixels in the following tests. 



2 4 6 8 10 

Overlap (pixels) 


Figure 11: Computation time with respect to the overlap. 
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There exist other methods without overlapping which may be theoretically faster 
m- In a forthcoming work, we will consider this class of methods. However, as 
we have seen, an overlap of 5 pixels is enough, so for large images the additional 
time can be neglected. On the figure we present the evolution of the estimation 
according to the iterations of the Schwarz method. We can see that since the third 
iteration, we already have a good estimation of the solution at the interfaces. 


Figure 12: Results for a 2x2 separation, one image per iteration of Schwarz and 
with an overlap of five pixels. 

In order to test the implementation of the multi-level parallelism, we have 
launched the same test case for different splittings of the image (2, 4 and 8 parts) 
and different numbers of CPU Per Part (CPP). In all cases, we have done ten iter¬ 
ations of the Schwarz method and we have kept an overlap of five pixels. On the 
figure we present the improvement of the computation times. 



2 3 4 5 6 7 8 

Number of parts 


Figure 13: Computation times of different configurations of the multi-level paral¬ 
lelism. 
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To understand the low efficiency of the multi-level implementation in the case 
where we use 32 CPU and split with 4 parts, we first represent on the figure the 
communication times obtained in the different cases. 



Number of parts 

Figure 14: Communication times with respect to the number of CPU and the 
splitting (does not include the intern communications due to the MUMPS solver). 

The increase of the communication times is not enough to explain the result 
obtained. Indeed, in every parallel implementations, the communication times in¬ 
crease with respect to the number of CPU used but it is usually balanced with the 
time saved in the computations. So, in order to give more details, we present in the 
figure \TE\ the times of the two main parts of the computation: the construction of 
the mass matrix and the LU factorization with the resolution part. We can see that 
adding more CPU to a part slightly increases the time of construction of the mass 
matrix. However, it allows to decrease the time of the LU factorization. On our 
architecture, if we split the images into several parts, the main part of the compu¬ 
tation is the mass matrix construction. If we use a large number of CPU per part, 
we reduce the difference between the time saved in the LU decomposition and the 
time lost in the matrix construction. More, we increase the communication times 
and the MUMPS parallel solver doesn’t balance that. 


14 

















Figure 15: Detailed computations times. 


To evaluate the accuracy of the computed flow in the case of varying illumina¬ 
tion, we use the RubberWhale test case where we have modified the first frame by 


increasing the brightness intensity (see figure 16). Since we haven’t modified the 
motion between the two pictures, the exact solution is the same that in the previous 
test case. 



Figure 16: Modified first frame of the RubberWhale sequence for the varying illu¬ 
mination test case. 

On the figure pTl we present the result obtained with and without the treatment 
of the brightness variation. On the right hand side, we can see that the solution 
without treatment is not well estimated. On the other hand, the model which uses 
the modified assumption allowing the brightness variation gives a much accurate 
optical flow still well oriented. 
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Figure 17: Left: Estimation with (right) and without (left) treatment of the varying 
illumination. The images were split in four parts. 


The next test (figure 18) consists of evaluating the efficiency of the adaptive 
regularization. The adaptation steps are done with the Schwartz iterations. We 
obtain a better definition of the edges. 


Figure 18: Estimated optical flow after ten iterations of adaptation (result for a 2x2 
separation). 

On the figure we can see that the adaptive regularization is still working 
for the test with non-constant brightness. In this case, the edges are even better 
defined. 



I 


Figure 19: Estimated optical flow after ten iterations of adaptation for the test with 
varying illumination (results for a 2x2 separation). 

Finally, on the figures and we have used two sequences of real images 
provided by the Centre d’etudes et d’expertise sur les risques, I’environnement, la 
mobilite et Vamenagement (Cerema). These pictures present two main interests. 
First, the high resolution. We have 2028 x 1098 pixels for the highway sequence 
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and 1524 X 1092 pixels for the wall sequence which corresponds to a vector field of 
about two billions of pixels to determine. The second interest is that the brightness 
and the texture are natural. 

The highway sequence presents a non-constant brightness with a more complex 
texture and a lot of occlusions as the white lines on the road or the large motion 
of the car. We can see on the figure the difference between the model with 
and without the varying illumination. On the left hand side, we have the solution 
without treatment. Again, we see that this is not correctly estimated because of 
the variation of the natural lighting. On the right, despite the occlusion areas 
which could be improved with a large displacement algorithm, we obtain a better 
approximation. 



Figure 20: Up: Test case of the Highway sequence given by the Cerema. Bottom: 
Optical flow estimation obtained with (right) and without (left) treatment of the 
varying illumination. The images were divided in sixteen parts. 

The figure represents a wall of a tunnel. On this sequence the motion to es¬ 
timate is linear. The challenge is to deal with irregular textures. This figure shows 
that our method is still efficient in this case. Indeed, the solution is quite smooth and 
linear except for the white square (bottom right) which may be improved by using a 
large displacement algorithm. This approach will be implemented in a further work. 

In the figure we give the computation times obtained with the two high 
definition tests. We recall that the pictures have about two billions of pixels and 
were not reduced. As in the previous cases, we perform ten iterations of the Schwarz 
algorithm. 
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Figure 21: Up: Test case of the Wall sequence given by the Cerema. Bottom: 
Optical flow estimation obtained with (right) and without (left) treatment of the 
varying illumination. The images were divided in sixteen parts. 



Number of parts 


Figure 22: Computation times of the high definition tests for different splittings. 


Conclusion 

In this article we have used the finite element method to solve the optical flow prob¬ 
lem with varying illumination. Thanks to the additive Schwarz method, we have 
implemented the domain decomposition in order to parallelize the computations. 
We have shown that this method is an efficient way to decrease the computation 
time and to handle high resolution sequences. We have also used a second level 
of parallelism in order to reduce the execution time even more and we have shown 
that such a massively parallel approach yields to an important gain of time. 
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The finite element method has also permitted to use an adaptive strategy of 
the regularization parameter, hence to efficiently combine the quality of the optical 
flow estimation with a method that preserves the edges and the fine features of the 
computed flow. 
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